car, ISLR, MASS, corrplotBoston{MASS}, Carseats{ISLR}lm{stats}, vif{car}연속형인 두 변수 사이의 종속성을 측정하는 도구인 상관계수(correlation coefficient)를 기초로 해 시행하는 분석 절차를 상관분석이라 함
두 변수 사이의 함수 관계를 규명하는 회귀분석의 전 단계에 산점도 분석과 함께 반드시 시행되는 절차임
상관계수는 Karl Pearson (1857 – 1936)이 제안한 개념으로 정의는 다음과 같음
\[ \rho = \frac{\textsf{Cov}(X, Y)}{\textsf{sd}(X)\textsf{sd}(Y)} \]
이 측도는 두 연속형 변수 사이의 선형 종속성(linear dependency)의 정도를 측정
\(-1\) ~ \(+1\) 사이의 값을 가지며, 0 근처이면 선형종속성이 유의하지 않음을, +1 또는 -1이면 완벽한 선형종속성을 가짐을 의미
스피어만(Spearman)의 순위상관계수는 변수의 순위 사이의 상관계수를 계산한 것으로 선형 종속성을 넘어 단조적(monotone) 종속성을 측정해 줌
켄달(Kendall)의 타우(\(\tau\))는 데이터에서 두 자료점을 선택해 concordant/discordant 여부를 따져서 concordant pair와 discordant pair 개수의 차이를 이용해 변수 간 연관성을 측정하는 도구임
round(cor(x1, x7, method = "spearman"), 4)
## [1] 0.7189
round(cor(x1, x7, method = "kendall"), 4)
## [1] 0.5167
round(cor(x10, x11, method = "spearman"), 4)
## [1] 0.5217
round(cor(x10, x11, method = "kendall"), 4)
## [1] 0.3453
iris 데이터pairs(iris[,1:4],
main = "Anderson's Iris Data -- 3 species", pch = 21,
bg = c("red", "green3", "blue")[unclass(iris$Species)])
setosa 품종 데이터만 추출해 변수들 간의 상관성을 분석Setosa <- subset(iris, Species == "setosa", select = -Species)
round(cor(Setosa), 4)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000 0.7425 0.2672 0.2781
## Sepal.Width 0.7425 1.0000 0.1777 0.2328
## Petal.Length 0.2672 0.1777 1.0000 0.3316
## Petal.Width 0.2781 0.2328 0.3316 1.0000
round(cor(Setosa, method = "spearman"), 4)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000 0.7553 0.2789 0.2995
## Sepal.Width 0.7553 1.0000 0.1799 0.2865
## Petal.Length 0.2789 0.1799 1.0000 0.2711
## Petal.Width 0.2995 0.2865 0.2711 1.0000
round(cor(Setosa, method = "kendall"), 4)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 1.0000 0.5973 0.2173 0.2311
## Sepal.Width 0.5973 1.0000 0.1426 0.2343
## Petal.Length 0.2173 0.1426 1.0000 0.2217
## Petal.Width 0.2311 0.2343 0.2217 1.0000
두 변수 사이의 상관관계가 있는 지 여부를 검정할 수도 있음.
with(Setosa, cor.test(Sepal.Width, Petal.Width))
##
## Pearson's product-moment correlation
##
## data: Sepal.Width and Petal.Width
## t = 1.6581, df = 48, p-value = 0.1038
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0487543 0.4800023
## sample estimates:
## cor
## 0.232752
with(Setosa, cor.test(Sepal.Width, Petal.Width,
method = "spearman", exact = FALSE))
##
## Spearman's rank correlation rho
##
## data: Sepal.Width and Petal.Width
## S = 14858, p-value = 0.04366
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.2865359
with(Setosa, cor.test(Sepal.Width, Petal.Width,
method = "kendall", exact = FALSE))
##
## Kendall's rank correlation tau
##
## data: Sepal.Width and Petal.Width
## z = 2.0593, p-value = 0.03946
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.2342674
회귀분석이란?
종속변수(dependent variable) 혹은 반응변수(response variable)
다른 변수에 의해 설명되는 변수
보통 측정오차를 수반함
독립변수(independent variable) 혹은 설명변수(explanatory variable)
다른 변수에 영향을 주거나 그 변수를 설명하는 변수
보통 측정오차가 없다고 가정
회귀모형 \[ y = f(x) + \epsilon \]
회귀분석의 활용
여러 통계적인 분석기법 가운데 가장 널리 쓰임
공학, 자연과학, 사회과학 및 의약학 등 모든 학문 분야에서 사용됨
분류 문제도 회귀분석 문제로 볼 수 있음
사례 : Advertising data
Questions
각 매체 광고비와 판매량간에 어떤 관계가 있는가?
관계가 있다면 얼마나 강해 보이는가?
광고비와 판매량 간의 관계가 선형인가?
반응변수 \(Y\)와 하나의 설명변수 \(X\)간에 선형함수관계를 가정하여 모형을 적합
실제 회귀함수가 선형이 아닐지라도 선형회귀는 해석 및 활용이 간편하다는 장점을 가짐
단순선형회귀모형
\[Y = \beta_0 + \beta_1 X + \varepsilon\]
기본가정 : 모형의 선형성, 오차항의 정규성, 독립성, 등분산성
최소제곱추정(least square estimate)
\[\mbox{RSS} = (y_1-\hat\beta_0 - \hat\beta_1 x_1)^2+(y_2-\hat\beta_0 - \hat\beta_1 x_2)^2+\cdots+(y_n-\hat\beta_0 - \hat\beta_1 x_n)^2\]
Measures of Fit: \(R^2\)
위 회귀모형에 따르면 \(Y\)의 변동의 요인은 설명변수 \(X\)와 오차 \(\epsilon\) 두 가지임
\(Y\)의 변동을 설명변수에 의한 변동과 오차에 의한 변동으로 분해할 수 있음 (분산분석, ANOVA)
\(R^2\)는 항상 0과 1 사이의 값을 가짐
\(R^2 = 0\)은 \(X\)를 이용해 설정했던 모형이 아무런 설명력이 없음을 의미하고, \(R^2 = 1\)은 그 모형으로 데이터를 완벽하게 설명한다는 것을 의미
set.seed(123)
N <- 50
x <- runif(N, min = -3, max = 3)
y <- 2 + 3*x + rnorm(N, sd = 5)
plot(x, y)
abline(2, 3, col = "navy", lwd = 3)
fit.lm <- lm(y ~ x)
plot(x, y)
abline(2, 3, col = "navy", lwd = 3)
abline(print(coef(fit.lm)), col = 2, lwd = 3)
## (Intercept) x
## 2.242598 3.318121
legend("topleft", c("True", "Fitted"), col = c("navy", "red"), lwd = c(3, 3))
summary(fit.lm)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.2789 -2.7893 -0.3284 2.7463 10.9306
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.2426 0.6665 3.365 0.00151 **
## x 3.3181 0.3804 8.723 1.82e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.702 on 48 degrees of freedom
## Multiple R-squared: 0.6132, Adjusted R-squared: 0.6051
## F-statistic: 76.08 on 1 and 48 DF, p-value: 1.822e-11
M <- 10000
beta1 <- rep(NA, M)
for ( r in 1:M ) {
y <- 2 + 0*x + rnorm(N, sd = 5)
beta1[r] <- coef(lm(y ~ x))[2]
}
hist(beta1, probability= TRUE, main = "", xlim = c(-3, 5),
col = "lightgray", border = "white", nclass = 100,
xlab = expression(beta[1]))
abline(v = coef(fit.lm)[2], col = 2, lty = 2)
points(coef(fit.lm)[2], 0, pch = 25, bg = "yellow", col = 2)
mean(beta1 >= coef(fit.lm)[2])
## [1] 0
adv <- read.csv(file = "./dat/Advertising.csv")
head(adv)
## X TV Radio Newspaper Sales
## 1 1 230.1 37.8 69.2 22.1
## 2 2 44.5 39.3 45.1 10.4
## 3 3 17.2 45.9 69.3 9.3
## 4 4 151.5 41.3 58.5 18.5
## 5 5 180.8 10.8 58.4 12.9
## 6 6 8.7 48.9 75.0 7.2
summary(adv)
## X TV Radio Newspaper
## Min. : 1.00 Min. : 0.70 Min. : 0.000 Min. : 0.30
## 1st Qu.: 50.75 1st Qu.: 74.38 1st Qu.: 9.975 1st Qu.: 12.75
## Median :100.50 Median :149.75 Median :22.900 Median : 25.75
## Mean :100.50 Mean :147.04 Mean :23.264 Mean : 30.55
## 3rd Qu.:150.25 3rd Qu.:218.82 3rd Qu.:36.525 3rd Qu.: 45.10
## Max. :200.00 Max. :296.40 Max. :49.600 Max. :114.00
## Sales
## Min. : 1.60
## 1st Qu.:10.38
## Median :12.90
## Mean :14.02
## 3rd Qu.:17.40
## Max. :27.00
par(mfrow = c(1, 3))
with(adv, plot(TV, Sales))
with(adv, plot(Radio, Sales))
with(adv, plot(Newspaper, Sales))
반응변수: Sales
독립변수: TV
회귀모형식: \[ \textsf{Sales} = \beta_0 + \beta_1 \cdot \textsf{TV} + \varepsilon \]
모형 적합
lm()을 사용하면 간단fit.lm <- lm(Sales ~ TV, data = adv)
with(adv, plot(TV, Sales))
abline(coef(fit.lm), col = "navy", lwd = 3)
분석 결과 정리
summary()를 이용해 요약summary(fit.lm)
##
## Call:
## lm(formula = Sales ~ TV, data = adv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.3860 -1.9545 -0.1913 2.0671 7.2124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.032594 0.457843 15.36 <2e-16 ***
## TV 0.047537 0.002691 17.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.259 on 198 degrees of freedom
## Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099
## F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16
설명변수가 TV 하나인 것보다는 Radio, Newspaper까지 다 사용하는 것이 어떨까?
다중회귀모형식 \[Y_i = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \cdots + \beta_{p}X_{p} + \varepsilon\]
\(\beta_0\)는 절편
\(\beta_j\)는 \(X_j\)가 한 단위 증가하고 나머지 \(X\) 변수의 값은 일정하게 유지될 때, \(Y\)의 평균 증가량에 해당 (partial regression coefficients)
최소제곱추정: 다음 식을 최소로 하는 \(\beta\) 계수
\[\text{MSE} = \frac{1}{n} \sum_{i=1}^n (Y_i-\hat{Y})^2 = \frac{1}{n} \sum_{i=1}^n (Y_i-\hat{\beta_{0}}-\hat{\beta_{1}}X_{1}-\cdots-\hat{\beta_{p}}X_{p})^2\]
fit.lm <- lm(Sales ~ TV + Radio + Newspaper, data = adv)
\(\beta_j\)가 0 인가?… 에 대한 답은 가설검정을 통해서 얻을 수 있음
\(\beta_j = 0\) 이면 \(X_j\)를 모형에 포함시킬 이유가 없음
각 계수별 \(t\)-test 결과 확인
summary(fit.lm)
##
## Call:
## lm(formula = Sales ~ TV + Radio + Newspaper, data = adv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8277 -0.8908 0.2418 1.1893 2.8292
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.938889 0.311908 9.422 <2e-16 ***
## TV 0.045765 0.001395 32.809 <2e-16 ***
## Radio 0.188530 0.008611 21.893 <2e-16 ***
## Newspaper -0.001037 0.005871 -0.177 0.86
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.686 on 196 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8956
## F-statistic: 570.3 on 3 and 196 DF, p-value: < 2.2e-16
\(X\) 변수들 중에서 적어도 하나는 유용한 예측 변수라고 확신할 수 있는가?
이 질문은 결국 \(\beta_1\) = \(\beta_2\) = \(\dots\) = \(\beta_p\) = 0 인가?… 와 같은 것임
\(F\)-test 결과 확인
모형 수정
Newspaper는 유의하지 않은 것으로 나타남
Newspaper를 제외한 나머지 두 개의 변수만을 포함한 회귀모형
fit.lm <- lm(Sales ~ TV + Radio, data = adv)
summary(fit.lm)
##
## Call:
## lm(formula = Sales ~ TV + Radio, data = adv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7977 -0.8752 0.2422 1.1708 2.8328
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.92110 0.29449 9.919 <2e-16 ***
## TV 0.04575 0.00139 32.909 <2e-16 ***
## Radio 0.18799 0.00804 23.382 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.681 on 197 degrees of freedom
## Multiple R-squared: 0.8972, Adjusted R-squared: 0.8962
## F-statistic: 859.6 on 2 and 197 DF, p-value: < 2.2e-16
예측 정확도 평가
Sales.pred <- predict(fit.lm)
plot(adv$Sales, Sales.pred, xlab = "Actual", ylab = "Predicted")
abline(0, 1, col = 2, lty = 2)
mean((adv$Sales - Sales.pred)^2) # mean(fit.lm2$residuals^2)
## [1] 2.78457
\[ \textsf{Sales} = 2.9211 + 0.0458 \cdot 100 + 0.1880 \cdot 45 \]
newx <- data.frame(TV = 100, Radio = 45)
predict(fit.lm, newx)
## 1
## 15.95632
확실한 방법은 all possible regression, 즉 모든 가능한 경우를 고려해 최적 모형을 찾는 것임
그러나 too many!!! (\(2^p\) with \(p\) predictors)
Forward selection (전진선택법) :
절편항만 포함하는 Null model로부터 출발하여 RSS 기준으로 변수를 하나씩 추가해 가는 방법
남은 변수가 모두 유의하지 않을 때까지 반복해서 시행
Backward elimination (후진소거법) :
모든 변수가 들어있는 Full model로부터 출발하여 가장 큰 p-value를 갖는 유의하지 않은 변수를 제거해 가는 방법
남은 변수가 모두 유의할 때까지 반복해서 시행
Stepwise regression (단계별 선택법) :
변수를 하나씩 추가하면서 새로 들어온 변수로 인해 유의하지 않은 변수는 제거해 가는 방법
새로운 변수의 추가가 없으면 중지
함수 step()을 이용하면 간단
step(fit.lm)
## Start: AIC=210.82
## Sales ~ TV + Radio
##
## Df Sum of Sq RSS AIC
## <none> 556.9 210.82
## - Radio 1 1545.6 2102.5 474.52
## - TV 1 3061.6 3618.5 583.10
##
## Call:
## lm(formula = Sales ~ TV + Radio, data = adv)
##
## Coefficients:
## (Intercept) TV Radio
## 2.92110 0.04575 0.18799
최적의 모델을 선택하기 위한 기준
ISLR 패키지의 Credit데이터를 사용library(ISLR)
summary(Credit)
## ID Income Limit Rating
## Min. : 1.0 Min. : 10.35 Min. : 855 Min. : 93.0
## 1st Qu.:100.8 1st Qu.: 21.01 1st Qu.: 3088 1st Qu.:247.2
## Median :200.5 Median : 33.12 Median : 4622 Median :344.0
## Mean :200.5 Mean : 45.22 Mean : 4736 Mean :354.9
## 3rd Qu.:300.2 3rd Qu.: 57.47 3rd Qu.: 5873 3rd Qu.:437.2
## Max. :400.0 Max. :186.63 Max. :13913 Max. :982.0
## Cards Age Education Gender Student
## Min. :1.000 Min. :23.00 Min. : 5.00 Male :193 No :360
## 1st Qu.:2.000 1st Qu.:41.75 1st Qu.:11.00 Female:207 Yes: 40
## Median :3.000 Median :56.00 Median :14.00
## Mean :2.958 Mean :55.67 Mean :13.45
## 3rd Qu.:4.000 3rd Qu.:70.00 3rd Qu.:16.00
## Max. :9.000 Max. :98.00 Max. :20.00
## Married Ethnicity Balance
## No :155 African American: 99 Min. : 0.00
## Yes:245 Asian :102 1st Qu.: 68.75
## Caucasian :199 Median : 459.50
## Mean : 520.01
## 3rd Qu.: 863.00
## Max. :1999.00
Gender와 같은 범주형변수를 어떻게 회귀식에 포함시킬 것인가?
Dummy variable (가변수)로 데이터 코딩
Gender와 같은 이항변수의 경우 Male이면 0, Female이면 1로 입력
Ethnicity 같은 다항변수의 경우에는 one-hot 인코딩: African American는 (1, 0), Asian은 (0, 1)
Balance를 반응변수로, Income과 Gender을 설명변수로 사용
이항변수인 Gender를 다음과 같이 코딩
\[\textsf{Gender} = \begin{cases} 0 & \quad \text{if male} \\ 1 & \quad \text{if female}\\ \end{cases}\]
\[\textsf{Balance} = \beta_{0} + \beta_{1} \cdot \textsf{Income} + \beta_{2} \cdot \textsf{Gender} = \begin{cases} \beta_{0} + \beta_{1} \cdot \textsf{Income} & \quad \text{if male} \\ \beta_{0} + \beta_{1} \cdot \textsf{Income} + \beta_{2} & \quad \text{if female}\\ \end{cases}\]
fit.lm <- lm(Balance ~ Income + Gender, data = Credit)
summary(fit.lm)
##
## Call:
## lm(formula = Balance ~ Income + Gender, data = Credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -791.23 -351.34 -51.57 328.18 1112.87
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 233.7663 39.5322 5.913 7.24e-09 ***
## Income 6.0521 0.5799 10.437 < 2e-16 ***
## GenderFemale 24.3108 40.8470 0.595 0.552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 408.2 on 397 degrees of freedom
## Multiple R-squared: 0.2157, Adjusted R-squared: 0.2117
## F-statistic: 54.58 on 2 and 397 DF, p-value: < 2.2e-16
Male) 대비 여성(Females)의 Balance 초과값의 평균을 의미\[x_{i1} = \begin{cases} 1 & \quad \text{if Asian} \\ 0 & \quad \text{if not Asian}\\ \end{cases}\]
\[x_{i2} = \begin{cases} 1 & \quad \text{if Caucasian} \\ 0 & \quad \text{if not Caucasian}\\ \end{cases}\]
\[Y_i = \beta_{0} + \beta_{1}\, x_{i1} + \beta_{2}\, x_{i2} + \varepsilon_i = \begin{cases} \beta_{0} + \beta_{1} + \varepsilon_i & \quad \text{if Asian} \\ \beta_{0} + \beta_{2} + \varepsilon_i & \quad \text{if Caucasian}\\ \beta_{0} + \varepsilon_i & \quad \text{if African American}\\ \end{cases}\]
\(X_1\)의 증가로 인한 \(Y\)의 (증가 혹은 감소와 같은) 효과가 \(X_2\)값의 수준에 따라 달라질 때, interaction이 존재한다고 함
예제: Advertising 데이터
TV와 Radio의 광고는 모두 판매량을 증가시킴
두 가지 미디어 광고에 돈을 사용하는 것이 한 가지 광고에만 사용하는 것보다 판매량을 더 많이 증가시킴. Why?
\[\textsf{Sales} = \beta_{0} + \beta_{1} \times \textsf{TV} + \beta_{2} \times \textsf{Radio} + \beta_{12} \times \textsf{TV} \times \textsf{Radio}\]
\[\textsf{Sales} = \beta_{0} + (\beta_{1} + \beta_{12} \times \textsf{Radio}) \times \textsf{TV} + \beta_{2} \times \textsf{Radio}\]
\[\textsf{Sales} = \beta_{0} + (\beta_{2} + \beta_{12} \times \textsf{TV}) \times \textsf{Radio} + \beta_{2} \times \textsf{TV}\]
fit.lm <- lm(Sales ~ TV * Radio, data = adv)
summary(fit.lm)
##
## Call:
## lm(formula = Sales ~ TV * Radio, data = adv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3366 -0.4028 0.1831 0.5948 1.5246
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.750e+00 2.479e-01 27.233 <2e-16 ***
## TV 1.910e-02 1.504e-03 12.699 <2e-16 ***
## Radio 2.886e-02 8.905e-03 3.241 0.0014 **
## TV:Radio 1.086e-03 5.242e-05 20.727 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9435 on 196 degrees of freedom
## Multiple R-squared: 0.9678, Adjusted R-squared: 0.9673
## F-statistic: 1963 on 3 and 196 DF, p-value: < 2.2e-16
TV에 돈을 더 사용하는 것은 평균적인 판매량을 0.0191 + 0.0011 \(\times\) Radio만큼의 빠르기로 증가시킴
Radio에 돈을 더 사용하는 것은 평균적인 판매량을 0.0289 + 0.0011 \(\times\) TV으로 증가시킴
library(plot3D)
x <- adv$TV
y <- adv$Radio
z <- adv$Sales
fit.lm1 <- lm(z ~ x + y)
fit.lm2 <- lm(z ~ x * y)
grid.lines <- 200
x.pred <- seq(min(x), max(x), length.out = grid.lines)
y.pred <- seq(min(y), max(y), length.out = grid.lines)
xy <- expand.grid(x = x.pred, y = y.pred)
z1.pred <- matrix(predict(fit.lm1, newdata = xy),
nrow = grid.lines, ncol = grid.lines)
fitpoints1 <- predict(fit.lm1)
z2.pred <- matrix(predict(fit.lm2, newdata = xy),
nrow = grid.lines, ncol = grid.lines)
fitpoints2 <- predict(fit.lm2)
scatter3D(x, y, z, pch = 18, cex = 2,
theta = 50, phi = 16, ticktype = "detailed",
xlab = "TV", ylab = "Radio", zlab = "Sales",
surf = list(x = x.pred, y = y.pred, z = z1.pred,
fit = fitpoints1, col = "white", shade = 0.05),
main = "Advertising: w/o interaction", col = "red", alpha = 0.5)
scatter3D(x, y, z, pch = 18, cex = 2,
theta = 50, phi = 16, ticktype = "detailed",
xlab = "TV", ylab = "Radio", zlab = "Sales",
surf = list(x = x.pred, y = y.pred, z = z2.pred,
fit = fitpoints2, col = "white", shade = 0.05),
main = "Advertising: with interaction", col = "red", alpha = 0.5)
Credit 데이터 \[
\textsf{Balance}
= \beta_{0} + \beta_{1} \times \textsf{Income} + \beta_{2} \times \textsf{Student}
\] vs. \[
\textsf{Balance}
= \beta_{0} + \beta_{1} \times \textsf{Income} + \beta_{2} \times \textsf{Student} + \beta_{12}\times \textsf{Income}\times\textsf{Student}
= \begin{cases} \beta_{0} + \beta_{1} \times \textsf{Income} & \quad \text{if}\,\textsf{Non-Student} \\
(\beta_{0} + \beta_{2}) + (\beta_{1} + \beta_{12})\times\textsf{Income} & \quad \text{if}\, \textsf{Student}\\ \end{cases}
\]fit.lm1 <- lm(Balance ~ Income + Student, data = Credit)
fit.lm2 <- lm(Balance ~ Income * Student, data = Credit)
summary(fit.lm1)
##
## Call:
## lm(formula = Balance ~ Income + Student, data = Credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -762.37 -331.38 -45.04 323.60 818.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 211.1430 32.4572 6.505 2.34e-10 ***
## Income 5.9843 0.5566 10.751 < 2e-16 ***
## StudentYes 382.6705 65.3108 5.859 9.78e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 391.8 on 397 degrees of freedom
## Multiple R-squared: 0.2775, Adjusted R-squared: 0.2738
## F-statistic: 76.22 on 2 and 397 DF, p-value: < 2.2e-16
summary(fit.lm2)
##
## Call:
## lm(formula = Balance ~ Income * Student, data = Credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -773.39 -325.70 -41.13 321.65 814.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 200.6232 33.6984 5.953 5.79e-09 ***
## Income 6.2182 0.5921 10.502 < 2e-16 ***
## StudentYes 476.6758 104.3512 4.568 6.59e-06 ***
## Income:StudentYes -1.9992 1.7313 -1.155 0.249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 391.6 on 396 degrees of freedom
## Multiple R-squared: 0.2799, Adjusted R-squared: 0.2744
## F-statistic: 51.3 on 3 and 396 DF, p-value: < 2.2e-16
with(Credit, plot(Income, Balance, type = "n", main = "Model without interaction terms"))
beta <- coef(fit.lm1)
abline(beta[1], beta[2], lwd = 3)
abline(beta[1] + beta[3], beta[2], lwd = 3, col = 2)
legend("topleft", c("Student", "Non-student"), col = c(2, 1), lwd = c(3, 3))
with(Credit, plot(Income, Balance, type = "n", main = "Model with interaction terms"))
beta <- coef(fit.lm2)
abline(beta[1], beta[2], lwd = 3)
abline(beta[1] + beta[3], beta[2] + beta[4], lwd = 3, col = 2)
legend("topleft", c("Student", "Non-student"), col = c(2, 1), lwd = c(3, 3))
다항회귀 polynomial regression
with(Auto, plot(horsepower, mpg, col = "gray", main = "Auto dataset"))
fit.lm1 <- lm(mpg ~ horsepower, data = Auto)
fit.lm2 <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
x <- seq(from = min(Auto$horsepower), to = max(Auto$horsepower), by = 0.0025)
with(Auto, plot(horsepower, mpg, col = "gray"))
pred1 <- predict(fit.lm1, newdata = data.frame(horsepower = x))
lines(x, pred1, col = "orange", lwd = 3)
pred2 <- predict(fit.lm2, newdata = data.frame(horsepower = x, x^2))
lines(x, pred2, col = "darkgreen", lwd = 3)
legend("topright", c("Linear", "Quadratic"), col = c("orange", "darkgreen"), lwd = c(3, 3))
Note: 비선형회귀 nonlinear regression
최소제곱추정법은 오차항의 등분산성을 가정할 때 최고의 방법임
그러나 종종 이 가정이 깨지는 경우가 발생
fit.lm <- lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
summary(fit.lm)
##
## Call:
## lm(formula = mpg ~ horsepower + I(horsepower^2), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.7135 -2.5943 -0.0859 2.2868 15.8961
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.9000997 1.8004268 31.60 <2e-16 ***
## horsepower -0.4661896 0.0311246 -14.98 <2e-16 ***
## I(horsepower^2) 0.0012305 0.0001221 10.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.374 on 389 degrees of freedom
## Multiple R-squared: 0.6876, Adjusted R-squared: 0.686
## F-statistic: 428 on 2 and 389 DF, p-value: < 2.2e-16
plot(predict(fit.lm), fit.lm$residuals, xlab = "Predicted", ylab = "Residuals")
abline(h = 0, col = 2, lty = 2)
다중공선성 문제 multicollinearity
강한 상관관계가 있는 설명변수들을 모형에 포함시키는 경우 최소제곱추정에 문제를 발생시킴
추정된 회귀계수값이 이상하거나 표준오차가 매우 커짐
극단적인 경우 모형 적합이 아예 이루지지 않을 수 있음
pairs(Auto[,-(7:9)])
fit.lm <- lm(mpg ~ displacement + horsepower + weight + acceleration, data = Auto)
summary(fit.lm)
##
## Call:
## lm(formula = mpg ~ displacement + horsepower + weight + acceleration,
## data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.378 -2.793 -0.333 2.193 16.256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.2511397 2.4560447 18.424 < 2e-16 ***
## displacement -0.0060009 0.0067093 -0.894 0.37166
## horsepower -0.0436077 0.0165735 -2.631 0.00885 **
## weight -0.0052805 0.0008109 -6.512 2.3e-10 ***
## acceleration -0.0231480 0.1256012 -0.184 0.85388
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.247 on 387 degrees of freedom
## Multiple R-squared: 0.707, Adjusted R-squared: 0.704
## F-statistic: 233.4 on 4 and 387 DF, p-value: < 2.2e-16
분산팽창계수 VIF
Variance inflation factor (VIF) \[ \textsf{VIF}(\beta_j) = \frac{1}{1 - R_{j | -j}^2}\] 단, \(R_{j| -j}^2\)는 \(x_j\)를 반응변수로 하고 나머지 변수들을 설명변수로 해서 적합시킨 선형회귀모형의 \(R^2\)
보통 VIF값이 10 이상인 설명변수가 있으면 다중공선성 문제가 있다고 봄
library(car)
## Loading required package: carData
vif(fit.lm)
## displacement horsepower weight acceleration
## 10.686922 8.823022 10.284456 2.603255
Non-normal responses: Logistic regression, Poisson regression?
Non-linearity: Kernel smoothing, splines and generalized additive models
Interactions: Tree-based methods, bagging, random forests and boosting (these also capture non-linearities)
Regularized fitting: Ridge regression and lasso
지금까지의 회귀모형은 예측 대상 변수인 \(y\)값에 대해 정규분포 가정을 할 수 있는 경우였다. 그러나 실제 분석에서는 예측 대상 변수가 범주형이거나 계수형이어서 정규분포 가정이 적절치 않은 경우가 태반이다. 이러한 경우를 위한 선형모형을 일반화선형모형(generalized linear model, GLM)이라 한다.
로지스틱 회귀모형은 반응변수 \(y\)가 성공(1), 실패(0) 두 가지 값(예: 환자 생존 여부, 전염병 발병 여부, 신용 부도 발생 여부, 고객 claim 발생 여부 등)을 갖는 이항(binary) 변수인 경우에 대한 회귀모형으로 일반화선형모형의 대표적 방법론이다.
반응변수가 연속형인 경우와 달리, 이항 변수인 경우 \(E(y|x) = P(y=1\,|\,x)\)가 되어 취할 수 있는 값의 범위가 0과 1 사이의 값으로 제한된다. 이를 무시하고 통상적인 회귀분석을 적용하면 예측값이 허용 범위인 0과 1 사이를 벗어나는 경우가 발생해 문제가 된다.
따라서 다음과 같이 \({\rm logit} (p) = \log\frac{p}{1-p}\)을 이용(이런 함수를 link function이라 함)해 변환한 모형을 고려한다.
\[
{\rm logit}P(y = 1\,|\,x_1,...,x_p) = \beta_0 + \beta_1 x_1 + ... + \beta_p x_p
\] 참고로 로짓함수의 역함수는 sigmoid 함수 \(\sigma(s) = \frac{1}{1 + e^{-s}}\)이다. 즉, \[
P(y = 1\,|\,x_1,...,x_p) = \sigma(\beta_0 + \beta_1 x_1 + ... + \beta_p x_p)
\] 으로 표현할 수 있다.
faraway 패키지의 built-in 데이터셋인 orings 데이터
1986년 우주왕복선 챌린저호 폭발 사건이 로켓엔진의 O-ring 불량과 관련이 있다고 알려져있다.
이 데이터셋은 이전 23회의 우주 왕복 임무 과정에서 수집된 데이터이다.
temp: 발사 때 기온 (화씨)damage: 6회 시험 중 O-ring 손상이 발생한 횟수
library(faraway)
##
## Attaching package: 'faraway'
## The following object is masked _by_ '.GlobalEnv':
##
## logit
## The following objects are masked from 'package:car':
##
## logit, vif
data(orings)
str(orings)
## 'data.frame': 23 obs. of 2 variables:
## $ temp : num 53 57 58 63 66 67 67 67 68 69 ...
## $ damage: num 5 1 1 1 0 0 0 0 0 0 ...
fit.glm <- glm(cbind(damage, 6 - damage) ~ temp, data = orings,
family = binomial(link = logit))
#fit.glm <- glm(cbind(damage, 6 - damage) ~ temp, data = orings,
# family = binomial(link = probit))
summary(fit.glm)
##
## Call:
## glm(formula = cbind(damage, 6 - damage) ~ temp, family = binomial(link = logit),
## data = orings)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9529 -0.7345 -0.4393 -0.2079 1.9565
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.66299 3.29626 3.538 0.000403 ***
## temp -0.21623 0.05318 -4.066 4.78e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 38.898 on 22 degrees of freedom
## Residual deviance: 16.912 on 21 degrees of freedom
## AIC: 33.675
##
## Number of Fisher Scoring iterations: 6
plot(damage/6 ~ temp, xlim = c(25, 85), ylim = c(0, 1),
col = "navy", pch = 20,
xlab = "Temperature", ylab = "Pr(damage)", data = orings)
x <- seq(from = 25, to = 85, by = 0.01)
beta <- fit.glm$coefficients
lines(x, sigmoid(beta[1] + beta[2]*x), col = 2, lwd = 3)
abline(h = 0.1, col = 2, lty = 2)
print(thre <- min(x[sigmoid(beta[1] + beta[2]*x) <= 0.1]))
## [1] 64.1
abline(v = thre, col = 2, lty = 2)
points(thre, 0, pch = 25, bg = "yellow", col = 2)
predict(fit.glm, newdata = data.frame(temp = 31), type = "response") # OMG!
## 1
## 0.9930342